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6. Three-dimensional results are presented for a clipped delta wing with leading-edge 
sweep of 50.5° with a circular-arc airfoil section and for an aspect ratio 5 
rectangular wing with a NACA 64A006 airfoil section. 

7. A conservative shock point operator was derived for use at a mesh point where the 
steady flow is supersonic while the flow at the next point downstream is subsonic. 
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2.0 INTRODUCTION 


The purpose of the work described in this report was to develop a means for calculating 
air forces for use in flutter analyses of three-dimensional lifting surfaces in the 
transonic flight regime. Flutter is not only a significant problem at transonic speeds, 
but it has also proved difficult to predict analytically. These difficulties result not only 
from the mathematical complexities of the equations but also from computer resources 
required by the repetitive nature of flutter analyses performed during vehicle design. 

Various methods are currently under study for predicting unsteady transonic air forces, 
ranging from the relatively expensive finite difference models including time 
integrations to economical approximate procedures based on linear theory. The 
procedure of this report is intended to be intermediate in terms of computer machine 
resource usage and is based on a finite difference method developed by Ehlers in 
reference 1. The assumption of small perturbations from a uniform stream near the 
speed of sound retains the necessary complexity for describing flows with local 
supersonic regions. The application of the perturbation velocity potential restricts the 
solution to weak shocks which, for thin wings of reasonably good design, is not too 
limiting an assumption. When the flow is steady, the resulting nonlinear differential 
equation reduces to the well-known transonic small perturbation equation studied by 
Murman, Cole, and Krupp (refs. 2, 3, and 4). The unsteady differential equation is 
simplified by considering the flow to consist of the sum of two separate potentials 
representing the steady and unsteady effects. The assumption of small amplitudes of 
harmonic oscillation leads to a linear differential equation for the unsteady potential 
with variable coefficients depending on the steady flow. The resulting air forces are 
thus superposable and may be directly used in conventional flutter analysis 
formulations. 

The effect of thickness is included in the steady flow analysis. The unsteady analysis is 
carried out for a wing of vanishing thickness but submerged in a velocity potential 
distribution resulting from the steady analysis. As formulated, the shock is assumed to 
be fixed by the steady flow. It is noted that shock motion could be included in a linear 
fashion by introducing the perturbations of the unsteady motion into the 
Rankine-Hugoniot relations. 

Generally, the results of applying this procedure, as reported in references 5 through 9, 
have been encouraging. First, correlation of finite difference solutions for flat plate 
configurations with corresponding results from linear theory has been good for both 
two- and three-dimensional configurations. For mixed flow, where the solutions for a 
NACA 64A006 airfoil were compared with experimental data from Tijdeman and 
Schippers (ref. 10), the pattern of the pressure distribution closely followed that 
observed experimentally; however, the analytical pressure levels were generally higher 
than the measured levels. The reason for the discrepancy between theory and 
experiment is not known, but the discrepancy may be due to boundary layer or 
separation effects, or both, or to unknown problems associated with the theory or with 
the pressure measurements. Thus, the correlation studies for the two-dimensional case 
have been inconclusive because of the lack of knowledge of viscous effects and, for the 
three-dimensional case, because of a lack of experimental pressure data. 
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A significant cause for concern in the practical application of this procedure has been 
stability problems with the relaxation procedures used to solve the sets of finite 
difference equations. These stability problems - which are a function of reduced 
frequency, Mach number, and the size of the finite difference region - severely limit the 
use of this method in flow regimes of most interest. Solution stability is thus a major 
topic of this report. 

Section 6.0 is devoted to a discussion of the accuracy of solutions from the finite 
difference model in comparison with subsonic solutions for the flat plate. In addition, 
results are presented for a low-aspect-ratio delta wing and a moderate-aspect-ratio 
rectangular wing. 

A parallel study using finite difference methods on the unsteady transonic flow problem 
has been conducted by Traci, Albano, and Farr (refs. 7, 8, and 9). The resulting 
procedure concentrates in a consistent manner on the low-frequency regime. Their 
derived equations do not include the cross product term consisting of the derivative of 
the unsteady velocity potential <p\ with respect to time and of the second derivative of 
the steady velocity potential <po with respect to the flow-wise coordinate. In most of 
their applications, the second derivative with respect to time is left out. However, the 
formulation of the finite difference equations, the handling of the boundary conditions, 
and the use of a column line relaxation solution procedure appear very similar to the 
procedure used here. 
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.norm stream near the 


3.0 SYMBOLS AND ABBREVIATIONS 


a Streamwise dimension of mesh region, also value of x at left hand side of 

one-dimensional interval 

b Root semichord of wing, also vertical dimension of mesh region, also 

value of x at right hand side of one-dimensional interval. 

E Maximum error quantity 

f(x,y,t) Instantaneous wing shape defined by z 0 = 8f(x,y,t) 

f(y) Function defining wing trailing edge 

fo Undisturbed wing or airfoil shape 

f] Unsteady contribution to wing or airfoil shape 

h Distance between mesh points in one-dimensional problem 

i, j,k x,y,z subscripts for points in the mesh 

i VT 

K Transonic parameter, (l-M 2 )/(M 2 e) 

lx, ly , Dimensions of element used in residual discussion 

M Freestream Mach number 

ORF Overrelaxation factor 

q ( d 2 !€ - i*>(y-l)<p 0xx 

R Wall porosity parameter, also vector length used in Appendix B 

U 0 Freestream velocity 

x 0 ,y 0 , z o Physical coordinates, made dimensionless with the root semichord. 

x, y,z Scaled coordinates (xq, ,uyo> f° r the three-dimensional problem; the 

scaled coordinates for the two-dimensional problem are x and y, with x 
again being the direction of fluid flow. 

x'l.y'j.z'! Variables of integration 

xi, xt e Coordinates of leading and trailing edges 
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Coordinate of wing tip 
Angle-of-attack 



Ratio of specific heats for air 

Jump in pressure coefficient 

Jump in </>j at plane of wing or vortex wake 

Jump in , at wing trailing edge 

Thickness ratio or measure of camber and angle of attack 

(8/M> 2/3 

wM/( 1-M~) 

Scale factor on y<> and z 0 , \x = S l// ’ 3 M 2//3 

Dummy scaled coordinates for two-dimensional problem 

Air Density 

Unsealed perturbation velocity potential 
Steady scaled perturbation velocity potential 
Unsteady scaled perturbation velocity potential 
Wake integral defined in equation (B-l) 

Angular reduced frequency (semichord times frequency in radians per 
second divided by the freestream velocity, cob/U) 
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4.0 FORMULATION AND SOLUTION 

A detailed mathematical derivation of the method for the solution of the unsteady 
velocity potential for the flow about a harmonically oscillating wing is presented in 
reference 1. The discussion here will be limited to a brief outline of the procedure for 
the two-dimensional flow. 

The complete nonlinear differential equation was simplified by assuming the flow to be 
a small perturbation from a uniform stream near the speed of sound. The resulting 
equation for unsteady flow is 

[K - (7 - l)<P t -(y + 1 ) <P X \ <P XX + ^yy " ( 2 ^xt + < PtO/ € = 0 (D 


where K = (1 - M 2 )/M 2 e, M is the freestream Mach number of velocity Uq in the 
x-direction, x and y are made dimensionless to the semichord b of the airfoil and the 
time t to the ratio b/Uo- With the airfoil shape as a function of time defined by the 
relation 

y 0 = Sf(x,t) 

the linearized boundary condition becomes 


<P y = f x (x,t) + f t (x,t) 


(2) 


The quantity 8 is associated with properties of the airfoil (such as maximum thickness 
ratio, camber, or maximum angle of attack) and is assumed small. The coordinate y is 
scaled to the dimensionless physical coordinate yo according to 

y = § 1/3 ]v|2/3yQ 


and e is given in terms of 8 by 

e = (5/M) 2 / 3 

The pressure coefficient is found from the relation 

C p = - 2 6(<P x+ * t ) 

The preceding differential equation is simplified by assuming harmonic motion and by 
assuming the velocity potential to be separable into a steady-state potential and a 
potential representing the unsteady effects. We write for a perturbation velocity 
potential 

(x,y) + (x,y)e* 0Jt (3) 


and for the body shape 


y 0 = Sf(x,t) = 5[f 0 (x) + f](x)e lcot 
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Since the steady-state terms must satisfy the boundary conditions and the differential 
equation in the absence of oscillations, we obtain 


[K- ( 7 + l)^o 1 ^0 = 0 

x u xx yy 


with 


y 


0 00» 


y = 0, -1 <x< 1 


(4) 


(5) 


On the assumption that the oscillations are small and products of <p\ may be neglected, 
equations (1) and (2) with the aid of equations (4) and (5) yield 


|[K-(t+ D^o x J *l x | 



(2icj/e) i/Jj + q </>! = 0 


(6) 


where 


q = co 2 /e - ico (7 - 1 ) i^o 


xx 


subject to the wing boundary conditions 

=fi (x) + icofi (x), y = 0, -1 <x< 1 (?) 

y x 1 

A computer program for solving the steady-state transonic flow about lifting airfoils 
based on equations (4) and (5) was developed by Krupp and Murman (refs. 3 and 4). The 
output of this program or a similar program can be used in computing the coefficients 
for the differential equation of the unsteady potential. The similarity of the unsteady 
differential equation to the steady-state equation suggests that the method of column 
relaxation used by Krupp for the nonlinear steady-state problem should be an effective 
way to solve equation (6) for the unsteady potential <pi. Note that equation (6) is of 
mixed type; being elliptic or hyperbolic whenever equation (4) is elliptic or hyperbolic. 
Central differencing was used at all points for the y derivative and all subsonic or 
elliptic points for the x derivatives. Backward (or upstream) differences were used for 
the x derivatives at all hyperbolic points. 


The boundary condition that the pressure be continuous across the wake from the 
trailing edge was found in terms of the jump in potential A <pi to be 


Ay>j - A<£]^ e -ico(x-Xf- e ) 


(9) 


where A^jte is the jump in the potential at x = xt e just downstream of the trailing edge 
and is determined to satisfy the Kutta condition that the jump in pressure vanish at the 
trailing edge. The quantity A<pj is also used in the difference formulation for the 
derivative ^i yy to satisfy continuity of normal flow across the trailing-edge wake. 
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For the set of difference equations to be determinate, the value of (p 1 or its derivative 
must be prescribed on the mesh boundary. Following Klunker (ref. 11), we found an 
asymptotic integral representation for the far-field <p 1 potential, and for the related 
pressure potential ^ix + i^i* Because of the difficulty with convergence of the integral 
over the wake for the integral equation of the velocity potential, upstream and 
downstream boundary conditions for the mesh were given in terms of the pressure 
potential <p x 4- ic o<pi, for which the wake integral can be integrated in closed form. The 
value of <pi was computed at one point on the upper boundary and at one point on the 
lower boundary; the points were conveniently chosen to facilitate rapid convergence of 
the wake integral. The values of ip± at other points on the upper and lower boundaries 
were found by numerically integrating the quantity <p + i axpi with respect to x. 

The numerical solution to the resulting large order set of difference equations may be 
obtained by using a relaxation procedure. The initial solutions were obtained by using a 
line relaxation procedure. Convergence is determined by monitoring ERROR, the 
maximum change in the velocity potential between iteration steps. ERROR is defined 
as the maximum value over all i and j of 

1J 

r 

where <p i i j (n) is the unsteady velocity potential for the nth iteration, <^i i j <n " 1) 
corresponding potential for the preceding iteration, and r is the relaxation factor. The 
solution was considered converged when ERROR lCf 5 . In some cases, particularly for 
finer meshes and for the pitch mode, convergence was considered complete when 
ERROR =sl0 -4 . 
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5.0 RELAXATION SOLUTION STABILITY 


As discussed in a preceding NASA report by the authors (ref. 5), significant stability 
problems were encountered with the relaxation procedures used to solve the finite 
difference equations. Generally, these procedures paralleled those successfully used for 
the steady-state problem. In essence, this meant sweeping through the mesh with a line 
relaxation procedure. When the line of points was parallel to the freestream, it was 
called row relaxation; when the line was perpendicular to the flow, it was called column 
relaxation. 

The characteristics of the solution instability are as follows: 

1. It occurs with the flow is purely subsonic as well as mixed with locally supersonic 
regions. Thus, the instability is not dependent on the presence or absence of 
transonic shock flow. 

2. It appears to be a function of Xj = o>M/(l - M 2 ) and the size of the finite difference 
area for the two-dimensional problem or volume for the three-dimensional problem. 
An analysis of the flat plate with a uniform mesh yields for the critical value of \ ly 
the value of Xj above which the relaxation solution is unstable, 

r 1 ii 

x i * * * i = 7r ~ + ^ (10) 

1 critical La- Kb 2 J 

where a is the streamwise dimension of the mesh region, b is the height, and K is 
the transonic parameter. 

3. For a given condition (say a fixed Mach number and finite difference point setup), 
as Xj was increased the rate of convergence decreased until the solution started to 
diverge. Thus, the actual value of Xj for which the solution first diverges is 
ill-defined, although it is generally in the neighborhood of the value given by 
equation (10). 

Some insight into the causes of the instability may be obtained by considering the 
Helmholtz equation into which the difference equation for the oscillating flow over a 
flat plate may be transformed, namely, 

X X X + Xyy + M 2 X=0 dD 

It is well known that solutions to the Helmholtz equation may not be unique for given 
types of boundary conditions on a closed region since eigenfunctions with real 
eigenvalues can occur; i.e., functions representing standing waves for which 
homogeneous boundary conditions occur on the boundary. For the rectangular mesh 
area of width b and length a, the first eigenvalue associated with solutions of the 
Helmholtz equation with Dirichlet boundary conditions is the critical value of Xj just 
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presented. In terms of the relaxation procedure, it was shown in reference 5 that 
solution of a relaxation problem of the form 

IA(X { )] {^} = {r} (12) 

converges only when [A] is positive definite, and this holds for the unsteady problem 
when k 1 is less than X lcr [tical* 

Integral equation solutions currently in use for the linearized subsonic unsteady 
solutions employ only the outgoing wave solution for the kernel function. Similarly the 
outgoing wave solution is used to define Klunker-type (ref. 11) boundary conditions on 
the outer boundary of the mesh region. Apparently the incoming wave solution is picked 
up during the numerical solution. 

Investigations to remove or moderate the relaxation solution stability problem may take 
any of several paths. The approaches discussed here include (1) modifying the boundary 
conditions with the hope that the numerical solution would pick up only the outgoing 
wave solution, (2) using a coordinate transformation so that boundary conditions in the 
physical plane at infinity could be applied to the outer boundaries of a finite mesh 
region, (3) replacing the iterative relaxation solution with a full direct solution and 
thus solving for all the unknown velocity potentials at one time, and (4) using an 
overlapping subregion concept. 

Approaches that have also been considered to some degree and have not proved 
successful include the following: 

• Artificial manipulation of elements in [A] in order to provide a better conditioned 
matrix. In particular, an attempt was made to shift the eigenvalues of [ A] by 
addition of a large diagonal matrix to [ A] . (Such an addition must, of course, be 
compensated for by appropriate modification of the right-hand side of the system.) 
This modification did not improve the stability or convergence properties of the 
solution method. Subsequent theoretical investigations revealed that such a 
modification is essentially equivalent to doing underrelaxation on the original 
system. 

• A sequential mesh refinement system based on the procedures discussed by Brandt 
in reference 12. 

• A mathematical technique for making [A] positive definite for values of ki above 

Xicritical by premultiplication by the conjugate transpose of [ A] . This procedure 

and some results are presented by Hafez, Rizk, and Murman in reference 13. Our 
experience has been essentially the same as they describe; (i.e., that the 
convergence rate in the relaxation solution of the [A * A] system is very slow and 
that a small value for the maximum difference between iterations does not imply 
that the last iteration is correspondingly close to the true solution. 
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5.1 VARIATIONS IN OUTER BOUNDARY CONDITIONS 


The Klunker-type boundary conditions defined <pi on the upper and lower boundaries 
and set <p ix + iaxpi on the upstream and downstream boundaries of the finite difference 
region. Since these boundary conditions apparently did not effectively sort out the 
incoming waves from the outgoing waves, alternative conditions were explored. These 
included using an outgoing radiation-type condition on all four boundaries and also a 
porous wall boundary condition on the upper and lower boundaries. The mathematical 
forms for these boundary conditions are summarized in table 1. The porous wall 
conditions could be varied to form either a "free jet” by making the porosity parameter, 
R, very large or a "solid wall” condition by making R small. In practice, the parameter 
is usually fixed by some empirical method for specific wind tunnel conditions but, for 
the current work, the interest is on how the stability of the relaxation solution may be 
dependent on its value. 

Table 1.— Equations for Boundary Conditions 


— 
Boundary conditions 

Boundary 

Equation 

1. Outgoing wave 

Upstream 

^1x- ,a V M 0 1 

= 0 


Downstream 

M 

*1x + lo, 1 

= 0 


Upper 

, . i^Vk" . 

= 0 


Lower 

Y 1 - u z 

= 0 

2. Porous wafl 

Upper 

*1x + iw *1 +^-0 1y 

= 0 


Lower 

* 1x + ic ^1 "j^ly 

= 0 


The pilot program was modified so that all six combinations of outer boundary 
conditions shown in table 2 could be run; that is, either of the two conditions on the 
upstream and downstream boundaries could be run with any one of the three boundary 
conditions specified for the upper and lower boundaries. The free-jet and solid wall 
boundary conditions also were programmed explicitly and thus could be applied without 
the need for fixing a value for R. The test example consisted of a two-dimensional airfoil 
of vanishing thickness oscillating in harmonic pitch at a Mach number of 0.9. For this 
case and for the mesh dimensions that are used, the reduced frequency above which 
relaxation solutions are expected to be unstable according to equation (10) is about 0.1. 
The examples were run for a very coarse mesh (17 x 10), and the overrelaxation factor 
(ORF) was varied to make sure the solution instabilities were not due to too large an 
ORF. 
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Table 2. — Types of Boundary Conditions 


Upstream and 

downstream 

boundaries 

Upper and lower 
boundaries 

1. Klunker 

1. Klunker 

2. Outgoing wave 

2. Porous wall 


Free jet (large R) 
Intermediate 
Solid wall (small R) 


3. Outgoing wave 


In summary, the results of the calculations showed that the alternate boundary 
conditions used did not significantly improve the convergence of the solution. In some 
cases, a slight increase in the value of reduced frequency was observed for which 
convergent solutions could be obtained. No combination of boundary conditions would 
provide solution convergence above a reduced frequency of 0.18. Since the exact values 
of co at which a relaxation solution stops converging and starts diverging cannot be 
exactly determined anyway, the results of this investigation were not considered 
promising. 


5.2 COORDINATE TRANSFORMATION 

A second concept explored in hopes of removing the relaxation solution stability 
problem was a coordinate transformation that permits the boundary conditions at 
infinity to be used on the boundaries of the finite difference region; that is, the physical 
region to infinity is mapped into the limited area of the finite difference mesh in the 
calculation plane. The particular form of transformation that was used is that suggested 
by Carlson (ref. 14) which, as he points out, allows for a physically realistic behavior of 
the solution at infinity. The physical plane is divided into three regions by 
perpendicular lines through the leading and trailing edges of the airfoil. (See fig.l.) The 
physical plane coordinates (x,y) are related to the calculation plane coordinates (£,tj) by 
the following relations: 

In region I where f=s-l 


x = -1 + tan£-y (£ + 1)J + tanjjy (£ + 1) 3 J 


In region II, where -1<£<1 
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In region III, where 


x = 1 + tan 


[f«-l)] + tan[f «-l)3] 


and 


, 7T 

y = tan r\ 

Two different boundary conditions were used. The first consisted simply of making <p— 0 
on all four boundaries; the second, of using the outgoing wave conditions discussed in 
the preceding section. Here, the outgoing wave condition was applied at the midpoint 
between the boundary and the point adjacent to the boundary. 

These changes did not solve the relaxation solution stability problem. For a given Mach 
number, for example, relatively little (if any) change was noted in values of reduced 
frequency at which the solution became unstable. 

It is of interest to note that the combination of the coordinate transformation and the 
outgoing wave boundary condition provided results for the flat plate which very closely 
matched corresponding data from the NASA subsonic air force program (refs. 15 and 
16). A comparison of results from using outgoing wave conditions together with the 
coordinate transformation is shown in figure 2. It should be noted that the former 
results are for a 42 x 30 mesh while the latter results are for a significantly coarser 28 
x 20 mesh. 


5.3 COMPLETE DIRECT SOLUTION 


A "semidirect” solution procedure was examined by the authors in reference 5. The form 
of the equation solved at that time was 


[A (Xj)] 




RC^H-’V 


where contained an element for each interior mesh point. In other words, there 

was still an iteration required to update the vector {R(<Pi <n-1) )} on the right-hand side. 
Although very efficient for the small meshes for which it was used (i.e., permitted by 
the in-core solution capability), it was subject to the same type of solution instability as 
the relaxation solutions. However, it is possible to rewrite the equation so that all 
unknowns are on the left-hand side and the solution may be calculated without 
iteration. Consider, for example, the two-dimensional problem for purely subsonic flow. 
The mesh is set up to have IMAX points in x-direction and KMAX in the cross-flow 
direction. The points are sequenced by column, upstream to downstream. The unknowns 
consist of the <pi s interior to the outer boundaries. Thus, the indice N of the point I,K is 


N = (KMAX-2) * (1-2) + K 


for 2=s I =5 IMAX-1 and 2 =5 K =5 KMAX-1 
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Mach number = 0.9 
Reduced Frequency, 0J= 0.06 
Pitch about leading edge 



Figure 2 .— Jump in Pressure Coefficient Across a Flat Plate 
Oscillating in Pitch, M = 0.9, oj = 0.06 


The general form of the equation following eq. (24) of reference 1 (with all terms moved 
to the left-hand side) for points adjacent to the boundaries is of the form 


D<Pi-ik + Ayjjk-i + B<pik+ C^ik-i + E<pj+ik - 0 

With the sequencing as indicated above, the five terms in the coefficient matrix are in 
the following column locations: 

n D = KMAX* (1-3) + K 
n A = KMAX* (1-2) + K - 1 
n B = KMAX* (1-2) + K 
n c = KMAX* (1-2) + K + 1 
n E = KMAX* (1-1) + K 
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The bandwidth of the matrix is equal to ng-nD+1 or 2* (JMAX-2) + 1. The 
boundaries present special problems. For example, the <pi s adjacent to the wake, in 
addition to the usual dependency, are also functions of eight values of <pi at mesh points 
in the vicinity of the wing trailing edge (see eq. (41), (42), (85), and (86) of ref. 1). This 
significantly increases the bandwidth of the coefficient matrix. The <pi s for points on 
the outer boundaries are, using Klunker-type boundary conditions, functions of the ^’s 
at all other interior points in the mesh if the volume integrals are retained (see eq (110) 
and (114) of ref. 1). If the volume integrals are not retained (and this is the usual 
procedure), the boundary ^j’s remain functions of the A<pi’s across the wing and wake. 
Use of the outgoing wave boundary condition limits the dependency of the <pi for any 
point on the outer boundary to the immediate vicinity of that point. The bandwidth of 
the coefficient matrix is thus determined by the number of points in the wake. This 
complete or full direct procedure should provide answers over the full range of values 
except for the specific values of X x for which the matrix [A(Xj)] is singular. 

This procedure was first tested with a one-dimensional problem. There was no difficulty 
in obtaining solutions near the singular points. Accuracy, however, as measured against 
the analytic answers, did present difficulties, which are discussed in detail in section 5. 

The full direct solution was also investigated for the two-dimensional problems. One 
formulation included the coordinate transformation and the outgoing wave boundary 
conditions discussed previously. Use of the latter significantly reduced the bandwidth of 
the [A] matrix over what it would have been had Klunker-type outer boundary 
conditions been used, thus increasing the number of mesh points that could be handled 
by the in-core solution routines. 

The resulting program was used on the sample problem of the airfoil of vanishing 
thickness oscillating in pitch. As with the one-dimensional program, no trouble was 
encountered in obtaining solutions at frequencies well above values that had proved 
critical for the relaxation solution. However, once the neighborhood of the critical value 
had been reached or exceeded, very poor correlation with corresponding solutions from 
the NASA subsonic unsteady flow program was obtained; that is, as the value Xj was 
increased from subcritical values to supercritical values, correlation with the NASA 
program went from very good to very poor. The characteristics of this lack of correlation 
are discussed in detail in section 6. 

The original direct solution package did not contain a pivoting capability. Since concern 
was expressed about numerical accuracy of the solution in the neighborhood of the 
matrix singularities, a solution routine including partial pivoting with equilibration 
was inserted in the program. Although it could be determined that pivoting was used 
during the solution, the results remained exactly the same to the number of significant 
digits retained. 

In summary, the full direct solution provides solutions at values of Xj above the critical 
value. The solutions do not correlate well with corresponding solutions from the NASA 
subsonic unsteady flow program and are thus not considered reliable. Although these 
solutions have been obtained using routines that include partial pivoting, the lack of 
correlation does not appear to be due to numerical problems inverting the matrices. The 
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problem may be due to the restriction to a relatively coarse grid because of a limitation 
of the in-core solution routine and/or to the type of boundary conditions. This seems to 
be borne out by the results from the study of the one-dimensional problem for which an 
error analysis is easy to obtain. This is discussed in detail in section 6. 

5.4 OVERLAPPING REGIONS 

As noted in the general description of the relaxation solution stability problem at the 
beginning of this section, the critical value of Xj is inversely proportional to the size of 
the finite difference region over which the solution is being calculated. This fact 
suggested the possibility of solving a sequence of small overlapping regions. The critical 
value of Xj for each subregion would be large, and perhaps continuity between 
subregions could be achieved by iteration. In a sense, the basic line relaxation 
procedure (whether accomplished by rows or columns) is the extreme limit of this 
concept with each column or row acting as a separate subregion. Experience has already 
shown us that this does not work. 

However, some examples have been run using column relaxation with the finite 
difference solution region divided into two and three subregions vertically and with 
several variations in the amount of overlapping of the subregions. All results have been 
discouraging and there appears to be little point in extending this investigation further. 
A typical example using three subregions will be discussed next. 

This example is for a flat plate (no mixed flow) and a relatively coarse mesh (17 x 10), 
but it should provide a good indication of how the concept of overlapping regions will 
work. The solution region was divided into three subregions in the streamwise 
direction. The location of the mesh points in the x-direction and the corresponding 
indices are shown in figure 3a. Figure 3b presents a sequence of convergence histories 
with the range of x-indices over which the column relaxation solutions were performed 
indicated below each pass. First, the solution for the complete region was calculated to 
check solution stability. As shown in figure 3b, with an ORF of 1.0, the solution at first 
converges slowly but after some 90 iterations has started to diverge. Using an ORF of 
1.7, the solution is quite unstable and shows a general divergence trend. Included in the 
same figure are the convergence histories for the subregions. The calculation is started 
off by converging the middle section, which converges very rapidly, an ORF of 1.7 was 
used for this and all succeeding calculations. Then, as shown, the other sections were 
converged in succession. The <pi distribution was saved after each subregion solution 
and used as a starting point for the next solution. The overall convergence of the 
system, as noted from the starting error for each subregion solution, is marginal at best 
and would require many more solution sequences to determine whether the overall 
trend is convergent or divergent. 

Finally, the pressure plots for three different stops along the solution path are shown in 
figure 3c. These distributions do not appear to be converging either. This example is not 
considered to be completely conclusive as to the worth of the overlapping subregions 
concept. It is typical of what we have experienced with other similar examples. We have 
found no evidence that this concept would provide a practical means to avoid problems 
arising from relaxation solution instabilities. 
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6.0 NUMERICAL ACCURACY FOR LARGER VALUES OF 

The accuracy of the finite difference procedure of this report may be discussed in several 
different contexts. Previous reports (by the authors in refs. 1, 5, and 6 and Traci et al., 
in refs. 7, 8, and 9) have included numerical examples, the results of which are 
compared either with the experimental data of Tijdeman and Schippers (ref. 10) or with 
other analytical data. These analytical data may be for strictly subsonic flow (flow at 
high Mach number over a flat plate) or for more detailed transonic calculations 
including full shock effects such as those by Magnus and Yoshihara (ref. 17). The 
discussion here concentrates on the relationship between the critical value of ki (critical 
in terms of relaxation solution stability) and the accuracy of the finite difference 
solutions relative to more exact linear solutions. The examples to be discussed do not 
include shock effects. 


6.1 THE ONE-DIMENSIONAL PROBLEM 

In order to gain insight into the unsteady transonic problem as formulated in this 
report, a one-dimensional version of the flat plate problem was investigated. The 
one-dimensional analog of the two-dimensional equation (6) for a flat plate may be 
obtained by dropping the tp i term. Dividing the resulting equation by K, we have 

^I xx ' 2i Xl M ^l x + X l 2 (1 ‘ M 2 ) i^i = 0 (1 ° ; 


where Xj = oiM/(l - M 2 ). 

The exact general solution to equation (13) is 


ipj(x) = C j 


iXi ( 1+M)x 
e 


+ C 2 


-iXi (l-M)x 
e 1 


(14) 


where Cj and C 2 may be determined once the boundary conditions (end conditions) are 
specified. The derivation of equations (13) and (14) along with a detailed discussion of 
the exact general solution is presented in appendix A. An approximate solution over an 
interval [a,b] may also be found by transforming equation (13) to a finite difference 
equation with the solution being obtained by either a full direct solution (similar to that 
discussed in sec. 5.3) or by a point relaxation procedure. 

The interest here is in comparing answers obtained from the finite difference solution 
with corresponding answers from the exact solution. For this, the maximum error 
quantity E for a given reduced frequency w K , is defined as 


E (co K ) = 


Max 

1=1, IMAX 


*1 T " *1 T 

1 exact finite difference 


(15) 
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The investigation is aimed at determining the effect of the kind of boundary conditions 
used on E(o>k). First, it is clear from the exact solution that the solution for a given 
reduced frequency w or is made up of components with two substantially different 
wavelengths. For a given finite difference mesh (a given number of mesh points and 
specified mesh spacing), it would be expected that the short wavelength component 
would be less accurately represented than the long wavelength component; that is, a 
solution made up predominately of the short wavelength component would be less 
accurately determined using a finite difference calculation than a solution made up 
predominately of the long wavelength component. This has indeed proved to be the case 
as shown by examples presented in figure 4. Here two combinations of Dirichlet and 
Cauchy boundary conditions were used to obtain solutions. The first was set up so that 
the solution would consist solely of the short wavelength component and is denoted by 
the A-symbols in figure 4; and the second, set up so that the solution would consist 
solely of the long wavelength component, is denoted by the O symbols. The error level 
for the long wavelength component is significantly lower than that for the short 
wavelength component. 

1.0 
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Figure 4.— Comparison of Error Curve for Long and Short 
Wavelength Solutions 
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Second, it is of interest to know how the error varies with frequency. An analysis of a 
similar equation was made by Fischer and Usmani in reference 18. The equation 
studied was of the form 


i// xx + Xl 2 i|/ = ° 


(16) 


and is simply related to our one-dimensional equation by the transformation 


iXi Mx 

<p\ = \pe 


(17) 


Application of their analysis, based on equally spaced mesh points, to equation (16) 
shows that for small values of hXj, where h is the distance between adjacent mesh 
points and Dirichlet end conditions 

Eh 2 X , 3 

K 1 sin [X] (b - a)] (18) 

for some constant E independent of the reduced frequency and mesh spacing. In view of 
the close relation between the <pi and ib equations, we would expect the error behavior 
in the finite difference solution to be similar in both cases. Equation (18) displays 
several interesting characteristics. For example, the predicted error is directly 
proportional to the square of the mesh point spacing h and the third power of or, for 
fixed Mach number, the third power of the reduced frequency o>. Also, the presence of 
the sin [Xj(b - a)] in the denominator of the equation introduces singularities in the 
error curve at values of <o (or Xj) for which Xi(b - a) = mr , n = 1,2,... These values of X* 
correspond to eigenvalues of the analytical solution (eq. (14)): i.e., are values of Xj for 
which there is no unique analytical solutiqn. Except near these singularities, the error 
curve as a function of Xj behaves like Xj times a slowly varying modulation factor. 
Thus over much of the range of Xj the error is essentially proportional to Xj 3 . Very near 
X j = 0, the error is of course essentially proportional to Xi 2 since x/sin x — 1 as x— 0. 
This region is of little interest to the eigenvalue analysis, however. In view of the close 
relation between <pi and ^ equations, it is expected that the error behavior would also 
be the same for tp\ . 

The presence of the singularities in the curve of ERROR versus reduced frequency is 
shown in figure 5 by the A-symbols. It would appear that the eigenvalues for the 
analytic system do not coincide exactly with the eigenvalues for the finite difference 
system, as noted by the distortions in the curves with which the points have been 
connected. The calculation was set up so that E(g>k) would be evaluated at five points 
between each analytic eigenvalue. The singular behavior is the result of the evaluation 
of Cj and C 2 from a set of simultaneous equations that are a function of the applied 
boundary conditions. This set of equations may be written in the form 

[aCAj)] {C} = { 7 } 
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Figure 5. —Comparison of Error Curves for Boundary Conditions 
Yielding All-Real and Complex Eigenvalues 

where a is a 2 x 2 matrix that is a function of \\ (or a>), C is the two-element column 
matrix made up of and C 2 . The forms of a and y are a function of the nature of the 
boundary (end) conditions; i.e., whether they are Dirichlet, Neumann, or Cauchy. 
Moreover, for certain values of Xj, the determinant of a will be equal to zero. These 
certain values are eigenvalues. For values of Xi that correspond to eigenvalues, the 
solution for Cj and C 2 is not unique; that is, for Xj equal to eigenvalues, there is no 
unique solution to equation (13). 

It is interesting to note that the values of X 1? which are eigenvalues of a, may be either 
all-real or complex depending on the nature of the boundary conditions. It is readily 
shown that Dirichlet conditions on both ends or Neumann conditions on both ends lead 
to all-real eigenvalues. However, for certain combinations, such as mixed conditions 
(Dirichlet on one end and Cauchy on the other), the eigenvalues may be made complex. 
Under these circumstances, we would not expect the violent peak and valley behavior of 
the error plots that result from the all-real eigenvalues. This is indeed confirmed with 
the results shown in figure 5 when the boundary conditions are such as to yield complex 
eigenvalues. 
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This problem was originally studied to see if it would shed light on the relaxation 
solution instability problem. In particular, it was of interest to see if relaxation 
solutions could be obtained for boundary (end) conditions for which the eigenvalues are 
complex. However, tests with a relaxation solution of the one-dimensional system have 
not converged and thus having complex eigenvalues does not seem to materially affect 
the convergence. 

In addition it was noted that equation (18) implied that the error was essentially 
proportional to Aq 3 or o> 3 . An example of this is shown in figure 6 where an error curve 
for an example in which the singularity behavior has been suppressed is compared with 
a curve proportional to o> 3 . The correlation between the two is very good. Also included 
is a curve that is proportional to co 4 , as predicted by a conventional truncation analysis 
of the finite difference equation. 



Reduced frequency, (x) 


Figure 6“. — Variation of Error Curve With Reduced Frequency 
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In summary, analysis and experiment of the one-dimensional equation show that the 
error from the finite difference solution is essentially proportional to n X^ , and thus 
the number of points has to be expanded (or more specifically the mesh spacing reduced) 
in proportion to the 3/2 power of the frequency in order to retain accuracy. The level of 
the error is determined by the boundary conditions and, in turn, determines the relative 
contributions to the solution by the long and short wavelength components. The 
relatively larger part the long wavelength component plays, the smaller the level of 
error. Superposed on this general error curve can be a series of peaks and valleys with 
the peaks centered around the values of Xj (or co) that correspond to the real 
eigenvalues of the system of finite difference equations and are dependent on the 
boundary conditions. If the eigenvalues are complex, the peak- valley behavior of error 
curve is suppressed. 

These results would indicate that, for certain choices of boundary conditions and 
sufficiently fine mesh spacing, adequately accurate results may be obtained in the 
two-dimensional case using a full direct solution method. 

6.2 TWO-DIMENSIONAL EXAMPLES 

As noted in section 5.3, a complete direct solution using outgoing wave boundary 
conditions permits obtaining solutions at large values of reduced frequency, and 
solution stability no longer is a problem. However, for the mesh sizes used, the 
correlation between the finite difference solutions and linear theory becomes very poor. 
Results are presented here for a two-dimensional airfoil of vanishing thickness 
oscillating in pitch in a freestream of M = 0.9. Under these conditions, relaxation 
solutions would be expected to be unstable at reduced frequencies (based on the 
semichord) above approximately 0.12 according to equation (10). Results were obtained 
using both the linear theory program and the finite difference program. Very good 
correlation between the two theories was obtained at co = 0.06 (see fig. 2), and very poor 
correlation at co = 0.3 as shown in figure 7. The correlation was significantly degraded 
even at co = 0.09 as shown in figure 8. To test whether this phenomenon was a function 
of Xi rather than co, the same problem was rerun at a Mach number of 0.4 with reduced 
frequencies so that the values of Xj were the same. Correlation between results from 
linear and finite difference calculations, as shown in figures 9 and 10, was good for 
c o = 0.6 (corresponding to co = 0.06 at M = 0.9) and poor at co = 0.9 (corresponding to 
co = 0.09). The results at co = 3, which are not shown, were very bad. Thus, the results 
from the full two-dimensional transonic problem (although with nonmixed flow) appear 
to follow the same pattern as the results from the very simplified one-dimensional 
example. Indeed, the poor results appear to be due to the same cause as the peaks in the 
error curve shown in figure 5, but this requires further study. In particular, since the 
true eigenvalues of this problem are not known, it is difficult to assess whether the 
higher frequencies tried are near eigenvalues without further investigation of the 
sensitivity of accuracy as a function of frequency. 

These results were checked using a direct solution routine incorporating partial 
pivoting with equilibration. The results were not changed, although it was possible to 
tell that the pivoting portion of the routine had been used. Thus, the errors encountered 
with the two-dimensional calculations do not seem to be due to numerical problems 
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resulting from ill-conditioned matrices. Increasing the number of mesh points in order 
to improve correlation was not feasible with available computer resources. Decreasing 
the number of mesh points would not have provided a realistic representation of the 
physical problem. 



Linear theory computation 

In-phase 

— — Out-of-phase 

Finite difference computation 

• — In-phase 

— • — Out-of-phase 


Figure 7.— Jump in Pressure Coefficient Across a Flat Plate 
Oscillating in Pitch , M = 0.9 , to = 0.3 
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Figure 8.— Jump in Pressure Coefficient Across a Fiat Plate 
Oscillating in Pitch, M = 0.9, <x> = 0.09 







Figure 10.— Jump in Pressure Coefficient Across a Flat Plate 
Oscillating in Pitch, M = 0.4, to = 0.89 
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7.0 THREE-DIMENSIONAL PROGRAM STUDIES AND 

ANALYSES 


Modifications to the three-dimensional program as described in reference 5 are 
described in section 7.1. The results of applying the resulting program to the NASA 
transonic unsteady pressure model of low-aspect ratio clipped delta planform are 
presented in section 7.2. Section 7.3 presents results for a moderate-aspect-ratio 
rectangular wing. 


7.1 THE THREE-DIMENSIONAL PROGRAM 

The pilot three-dimensional program as described in reference 5 was restricted to lifting 
surfaces with rectangular planforms. This program has been revised, and its design and 
usage documented in reference 19. The revised program is valid for wings with 
aft-swept leading and trailing edges. The leading edge may be curved (of arbitrary 
shape), but the trailing edge must be straight. This last limitation is due to the method 
of programming rather than being a restriction on the theory. In addition, the program 
has been revised to: 

1 . Include the capability for row relaxation as well as the original column relaxation 

2. Make use of the anti-symmetry characteristics of the unsteady flow about 
symmetric wings so that only half the flow is actually calculated 

Row relaxation proved much faster than column relaxation for the two-dimensional 
problem. The same appears to be true from the minimal number of three-dimensinoal 
examples we have run. However, it should also be noted that solution instabilities have 
again been encountered in the mixed flow case, and the results of the following section 
for the configurations with thickness were obtained using column relaxation. It was 
noted in reference 5 that, for the two-dimensional problem, row relaxation was much 
more efficient than column relaxation in terms of reaching a specified degree of 
convergence in a minimum number of iterations. It was determined that in using row 
relaxation for mixed flow, additional terms must be included in the finite difference 
equation for hyperbolic points to avoid solution instabilities. These additional terms 
have not proved enough to avoid instabilities in the three-dimensional row relaxation 
solution, and it is assumed that the two-dimensional analysis of reference 5 should be 
extended to the three-dimensional equations. 


A derivation of the wake integral for a straight trailing edge perpendicular to the wing 
root was given in appendix B of reference 5. A general form, valid for wings with 
trailing edges, that may be described by a single valued function of the form 


x t (y) = 1 + f (y) 
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where f(y) ^=0, is derived in appendix B of this report. The resulting form again makes 
use of Gau8S-Laguerre integration and is directly parallel to the form derived in 
reference 5 for the trailing edge of an unswept rectangular wing. 

The three-dimensional program, which is considered to be a pilot program, has been 
provided to NASA-Langley and is documented in reference 19. The program has been 
modified and permits calculations including swept leading edges while using an 
unswept rectangular mesh point array. 

7.2 RESULTS FOR A DELTA WING 

This section presents the results of applying the pilot three-dimensional program to a 
wind tunnel model built by NASA-Langley for testing in the Langley Transonic 
Dynamics Tunnel. The model has a clipped delta planform with a 50.5 swept leading 
edge and a circular-arc profile with a thickness ratio of 6%. The model geometry is 
shown in figure 11. The model is designed to be oscillated in pitch and flapping, and 
every effort has been made to minimize the structural deflections resulting from these 
rigid body motions. The model is half-span and is mounted on the side of the tunnel 
through a splitter plate designed to remove the wall boundary layer. 

The calculations were performed at M = 0.9 for the wing oscillating in pitch and 
flapping at a reduced frequency based on the root semichord of 0.06. 

The steady-state pressure distribution for the wing is shown in figure 12. It was 
calculated using a program developed at NASA-Ames by Ballhaus and Bailey (ref. 20) 
and modified by The Boeing Company. It does not include a shock point operator. The 
ideas of Schmidt (ref. 21) were used to set up the mesh along the swept leading edge. 
The calculations were made for a mesh with 55 points in the flow direction, 32 points in 
the spanwise direction (half-span), and 36 points in the vertical direction. Convergence 
for the pitch mode and the flat plate configuration with ERROR ^ 10~ 4 and using row 
relaxation was obtained in about 100 iterations. Starting with this solution and using 
column relaxation, about 50 iterations were needed to obtain the solution for the 
circular-arc airfoil shape. With solutions calculated assuming symmetry with respect to 
the x-y plane and using a CDC 6600 computer with an FTN compiler, the number of 
CPU seconds per iteration was about 7 and the number per far-field update was about 9. 

The jump in pressure coefficient due to harmonic pitch and flapping is presented in 
figures 13 and 14. In each case, three different results are presented. The first result is 
from using the NASA subsonic unsteady three-dimensional airloads program (refs. 15 
and 16). This should compare directly with the second set of results, which are the finite 
difference results for a flat plate. The third data set is from using the finite difference 
program for the wing with the coefficients of the differential equation obtained from the 
nonlinear steady-state solution from the transonic small perturbation theory. 

Generally, linear results correlated very well with the corresponding finite difference 
results for a flat plate. This was particularly true for the pitch mode and only slightly 
less so for the flapping mode. Note that the scale used for the real part of the flapping 
mode is significantly larger than the scale used for the imaginary part. The failure of 
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the finite difference solution to provide a singularity at the wing leading edge at the 
root is attributed to the relative sparsity of points over the apex of the wing. In setting 
up the finite difference pattern, it was decided to emphasize the points on the aft 
portion of the wing and in the wake. From practical considerations, then, the planform 
apex was somewhat slighted in terms of points. 

The finite difference solution with thickness showed the usual peaks in the unsteady 
pressure in the region of the shocks as calculated in the steady flow. No experimental 
data are available at this time for comparison purposes. 

Results for the delta wing at an angle of attack are presented in figure 15. The first set 
of results was obtained using the nonlinear steady-state finite difference program for an 
angle of attack of 1.5°. The results are presented as jump in pressure coefficient per 
unit radian. The results from the nonlinear steady-state program were not fully 
converged; however, estimated converged results were indicated in the neighborhood of 
the shock as obtained using the Aitken-Shanks nonlinear transformation (o -process) of 
reference 22. The second set of results was from the unsteady program using a pitch 
mode and a very small reduced frequency of co = 0.00001. Only the real part of the 
resulting presure vector is plotted. The thickness effects in the unsteady program 
resulted from the steady velocity potential at zero angle of attack from the nonlinear 
program. The computer resources required to obtain the set of results from the linear 
unsteady program were significantly less than those required for the results from the 
nonlinear steady program. 

The pressure coefficient distributions from the two solutions exhibited the same 
characteristics with greater amplitude in the shock region for the linear unsteady 
solution for a very small frequency than for the steady nonlinear solution. Both 
solutions were obtained without using a shock point operator. 

7.3 RESULTS FOR A RECTANGULAR WING 

The revised three-dimensional program was also used to recalculate the pressure 
distribution over an aspect ratio 5 rectangular wing oscillating in harmonic pitch. A 
Mach number of 0.875 was used with a reduced frequency based on the root semichord 
of 0.06. These results as presented in references 5 and 6 were calculated using an 
incorrect scale factor on the steady-state velocity potential distribution. The effect of 
correcting this scale factor is to provide a noticeably larger pressure rise due to the 
presence of the shock. There still remains, however, a significant attenuation of this 
rise in going from the two-dimensional to the three-dimensional configuration. 

A mesh of 44 points in the flow direction, 32 points in the spanwise direction for the full 
span, and 26 points in the vertical direction was used. The finite difference region 
extended about one chord length in front of the leading edge and behind the trailing 
edge, about seven chord lengths above and below the wing surface, and slightly more 
than a semispan beyond the wingtip. The rerun has permitted a comparison of running 
times between the original program, using a KRONOS 2.1 operating system on the CDC 
6600 using the RUN compiler and the current program using an FTN compiler. The 
average number of CPU seconds per iteration is now approximately 2 compared to about 
8 before, and approximately 2V6 CPU seconds per far-field update compared to 9 before. 
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For the case shown, the converged solution (in this case the ERROR of eq. (9) was to be 
less than 10 4 ) was of the order of 180 iterations with the initial unsteady velocity 
potential distribution set to zeros. 

Figure 16 shows the steady-state pressure distribution for a NACA 64A006 profile, 
which was obtained using a program developed by Ballhaus and Bailey (ref. 20). The 
jump in pressure coefficient due to harmonic pitch about the planform leading edge is 
shown in figure 17, with three different results presented. The first results are from the 
NASA subsonic unsteady three-dimensional airloads program using linear theory (refs. 
15 and 16). These should compare directly with the second set of results calculated 
using the finite difference program and a flat plate airfoil section. The third data set is 
from using the finite difference program with the steady velocity potential distribution 
from the nonlinear steady-state solution for the wing with a NACA 64A006 profile. In 
addition, a two-dimensional result from finite difference theory for the same airfoil 
section is shown in the planform root plane. 

Generally, the linear results correlate very well with the corresponding finite difference 
results for a flat plate. The results including thickness display the pressure rise in the 
neighborhood of the shock that has been characteristic of corresponding experimental 
measurements (for example, see ref. 10). The three-dimensional results show a 
significant softening of the pressure rise in comparison with the two-dimensional 
results. Of concern is the apparent intensifying of the shock effect at the midpoint of the 
semispan of the wing. The reason for this result, which is not expected physically, is 
currently attributed to the way the finite difference operators are handled. The program 
is written to use central differencing for subsonic points (as determined from steady 
flow) and backward differencing for supersonic points. An abrupt change in the pattern 
of subsonic and supersonic points occurs on the chord adjacent to the one with the 
sharpest shock effects. 

In an attempt to smooth out the shock effects spanwise, a shock point operator in 
conservation form was introduced into the three-dimensional program. The derivation of 
the operator is given in appendix C. The result of using the shock point operator was to 
(1) significantly increase the effect of the shock on the unsteady pressure distribution 
and (2) smooth out the spanwise pressure distribution in the neighborhood of the shock. 
A comparison of distributions calculated with and without the shock point operator is 
shown in figure 18. Note the significant increase in the magnitude of the pressure rise 
due to the shock at the wing root, with a much smaller increment in the rise at 
midspan. No experimental data are available at this time for comparison purposes. 
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8.0 SUPERSONIC FREESTREAM 


Of significant interest is the inclusion of transonic flow effects in the calculation of 
oscillating air forces where the freestream flow is slightly supersonic. Of particular 
interest to the current work is whether or not the relaxation solutions become unstable 
in the same fashion when the freestream is supersonic as when it is subsonic. 


The differential and finite difference equations are the same for both the subsonic and 
supersonic freestream cases. The flow characteristics are sketched in figure 19, which 
shows the boundary conditions that were used in a pilot two-dimensional program. The 
unsteady velocity potential at the upstream boundary is set to zero. Since the flow is 
supersonic at the downstream boundary and backward differencing is used in the 
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Figure 19.— Boundary Conditions for Problem With Supersonic Freestream 



supersonic regions, boundary conditions need not be specified at the downstream 
boundary. Porous wall boundary conditions were convenient to use on the upper and 
lower boundaries. In practice, however, these boundaries should be set far enough out so 
that they do not affect the flow over the wing, and thus the pressure is independent of 
the porosity factor. 

As discussed by Traci et al. (ref. 9), the flat plate problem in which the steady-state 
velocity potential is constant may be solved by a single downstream pass with the 
relaxation procedure since nowhere in the flow is any point affected by points in the 
downstream columns. The problem of mixed flow with the pocket of subsonic flow buried 
within the supersonic flow is quite a different matter. Traci et al. noted relaxation 
solution instabilities in the neighborhood of M — 1.0 and, for the supersonic case, 
obtained two-dimensional solutions at M = 1.10 but not at M = 1.05. A priori, one may 
suspect that the finite subsonic region will have properties similar to the finite mesh of 
the subsonic freestream case, which results in instabilities in the relaxation process. 

In practice, numerical examples do not appear to admit such a simple explanation. A 
circular-arc airfoil was analyzed at two Mach numbers, M = 1.05 and 1.15. A simple 
pitching oscillation was studied. Some of our results have the characteristics of 
converging for a number of iterations and then diverging. Here the maximum difference 
between <p\ for successive iterations was used as a measure of convergence. If the 
convergence criteria were met before the divergence started, one would assume that one 
had obtained a valid solution. Under these circumstances, the use of overrelaxation 
factors (ORF) and underrelaxation factors (URF) other than unity increased the 
tendency for divergence. Hence, the calculations were run with ORF = URF - 1.0. The 
net result was that M = 1.15, with a relatively small subsonic region, the convergence 
characteristics were improved by raising the reduced frequency. At M = 1.05 with the 
attendant large subsonic region about the airfoil leading edge, convergence was 
improved by decreasing the reduced frequency. This latter behavior is what would be 
expected from experience with the subsonic freestream problem. 

These examples were rerun using the shock point operator of appendix C. Use of the 
operator noticeably improved the convergence characteristics at both Mach numbers but 
did not eliminate the relaxation solution instabilities. 

A summary of convergence experience with the supersonic freestream is given in table 
3. The table includes runs both with and without the shock point operator and includes 
a general description of where the maximum ERROR occurred for both converging and 
diverging examples. Since these calculations have been made with a limited number of 
variations in parameters such as the location of farfield boundaries, the number and 
spacing of mesh points, and the location of mesh points with respect to the sonic lines 
and subsonic regions, it is felt that firm conclusions are, as yet, unwarranted. 

There appear to be stability problems with the relaxation process in the supersonic 
freestream problem as well as with the subsonic problem. We suspect both have the 
same origins; that is, the eigen characteristics of the problem. However, numerical 
examples with the supersonic freestream problem do not give consistent convergence 
divergence behavior at M = 1.05 and M = 1.15. It is assumed that a full direct solution 
as described previously would provide solutions, but this has not been tried. 
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Table 3.— Comparison of Convergence Characteristics of Supersonic Free-stream 
With and Without Shock Point Operator ( SPO ) 


Reduced 

frequency, 

CO 

M = 1.05 

M = 1.15 

With SPO 

Without SPO 

With SPO 

Without SPO 

0.015 

Converges 

♦ Subsonic area 

* 3.5E-4 

Converges 

♦ In front of subsonic 
region off wing 

— 

— 

0.03 

Converges 
* Aft sonic line 
^3.4E-4 

Diverges 

♦ Front sonic line 

Diverges 

* Peak of subsonic 
region 

Diverges 

* Peak of subsonic 
region 

0.06 

Diverges 
* Aft sonic line 

Diverges 
♦ Aft sonic line 

Converges 

♦ Aft of subsonic 
region on wing 

* 2.0E-4 

Diverges 

♦ In front of subsonic 
region off wing 

0.12 

— 

— 

Converges 

♦ Aft of subsonic 
region on wing 

* 1.6E-4 

Diverges after 
35 iterations 
♦ In front of subsonic 
region off wing 

0.24 

— 

— 

Converges 

* Aft of subsonic 
region on wing 

* 0.8E-4 

Diverges after 
40 iterations 
♦ Front sonic line 


ORF = URF = 1.0 

& Maximum error after 100 iterations 
♦ Location of maximum error 
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9.0 PROPERTIES OF THE RESIDUAL 


It is customary for both steady and unsteady transonic finite difference solutions to use 
ERROR rather than the residual as a criterion for relaxation solution convergence. 
ERROR is the maximum difference between successive potential distributions divided 
by the relaxation factor as defined in equation (9). The residual is the difference 
between the right-hand side and the system matrix times the current approximation 
when the set of equations is written A<pj=R. This section summarizes a brief study of 
the characteristics of the residual with respect to the unsteady problem. 

First, it is noted that numerical examples, which are presented in appendix D show 
that a residual of the same order of magnitude as ERROR may be obtained by scaling 
the residual value with an associated area. For the purposes of this report, the 
RESIDUAL at a point ij is defined as the product of the value of the difference equation 
by a term proportional to the local mesh area, viz. 

( x i+l -x i _ 1 )(y j+1 -yj_i) 

4 

It is shown in appendix D that RESIDUAL may be interpreted as the excess (or deficit) 
of the mass flux within each mesh in the flow field. Since for an exact solution of the 
difference equation, this should be zero, it is also a measure of how close the relaxation 
solution is to being converged. 

It can be shown that RESIDUAL and ERROR curves, when plotted on an 
iteration history curve, should be essentially parallel to each other. A mathematical 
explanation of this phenomenon is also presented in Appendix D along with numerical 
examples for illustration purposes. 

Generally, it is not convenient to use the residual (or, for that matter, RESIDUAL) as a 
convergence test since evaluation requires a separate pass after all the velocity 
potentials have been updated. 
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10.0 AEROELASTIC ANALYSES 


Calculating air forces for flutter analyses can become very expensive when using the 
more complex aerodynamic procedures such as the method of this report. This is due to 
forces being functions not only of Mach number (as for steady-state analyses) but also of 
reduced frequency and due to the need to calculate pressures for a set of modes 
(generalized coordinates), which may number from 10 to 20 or more for low-aspect-ratio 
configurations. Also, flutter analyses are often required for small perturbations in mass 
and/or stiffness distributions from the basic configuration and, in general, this means 
recalculation of the generalized air force matrix. The question arises as to how the 
recalculations of the air force matrix can be handled most efficiently with respect to the 
procedures of this report for unsteady transonic flow. 

It is first noted that the basic differential equations are linear with spatial varying 
coefficients. The resulting air forces are thus superposable and may be directly used in 
conventional flutter analysis formulation. 

Next, the two kinds of numerical solutions to the finite difference equations that have 
been discussed are the line relaxation procedure and full direct solution. The former, 
which permits the solution to be calculated in sequences, is preferred because of the 
large number of finite difference points (and thus the large number of equations) even 
for two-dimensional problems. Indeed, it may well represent the only practical solution 
method for three-dimensional analyses. However, line relaxation does have instability 
problems for larger values of It was concluded in section 5 that for combinations of 
Mach number and reduced frequency where relaxation solutions are unstable, the most 
feasible method is the full direct solution. In matrix form, both of these procedures are 
written as 

mxm mxl mxl 

[A (M, co)] } = {R} l20/ 

It is assumed that outgoing wave boundary conditions are used on the far-field 
boundaries so that {R} does not depend on <p{ s. Also, the matrix sizes have been 
indicated above the equation. Here, m is the total number of finite difference points 
interior to the outer boundaries. The matrix {R} is a function of surface deformation 
(the mode shape) so that equation (20) may be rewritten as 

mxm mxl mxn nxl (21) 

[A (M, cu)]{^ 1 ( s )} = [T DW (u;)](f< s )} 


The matrix [Tdw(o>)] calculates the boundary conditions on the right-hand side from 

the modal matrix {f/ S) }. Here, n is the number of aerodynamic control points on the 
airfoil or wing and the superscript s denotes the mode shape. The size of n is expected to 
be on the order of 40 for the two-dimensional problem and on the order of 300 for the 
three-dimensional problem. 
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With respect to the direct solution, it is easy to conceive of an influence coefficient 
matrix that would be independent of structural characteristics; that is, a matrix which 
when postmultiplied by a modal matrix and premultiplied by its transpose would 
provide a set of generalized air forces directly. For example, equation (21) may be 
rewritten as 


mxl mxm mxn nxl 

(*>l (s) } - [A (M, co]" 1 [T dw (co)] j ^ ] 

The pressure is obtained by operating on the <pi’a, 

nx 1 nxm mxm mxn nx 1 

{ Ap(s>} = [Tp] [A (M, co] -1 [T dw (co)] |f< s > } 

The transformation matrix [T p ] transforms the velocity potential to the pressure 
distribution, and the generalized air force Q rs is found from 

lxn nxn / nxl » 

Qrs = Lf <r) J [T^Ap^)) (24) 

where the matrix [Tj ] performs the necessary integration of the product of the pressure 
due to mode s times the r mode shape. Substituting equation (23) in equation (24), we 
have 


( 22 ) 


(23) 


Qrs 


lxn nxn nxl 

l/ r) J([Tj] [T p ] [A (M, co)] -1 [T DW (co)]) ( f< s ) j 


(25) 


and the matrix product enclosed by the parentheses is just what is desired and is a 
complex matrix of order nxn. This represents a very manageable matrix in terms of size 
and number of operations to obtain Q rs . 

However, the critical problem is the size and banding characteristics of the matrix 
[A(M,w)]. For example, for practical two-dimensional problems, m will be 
approximately 1000 to 2000 and for practical three-dimensional problems it will be 
approximately 25 000 to 50 000. Also, [A(M,o>)] is complex, which essentially doubles 
its storage requirements. Here significant advantage can often be taken of the special 
case of flow symmetry with respect to the plane of the wing. While this may result in a 
two-dimensional problem of manageable size, it does not appear to do the same for the 
three-dimensional problem. An important characteristic of [A(M,o>)] is that it is banded; 
therefore, both [A(M,oO] and its LU decomposition form can be stored in significantly 
less space than the complete matrix. Thus, in practice, the inverse of [A], which would 
be a dense matrix, is not calculated. The unknowns that are found are the pi's of 
equation (21) or (22), and the solution is found from the LU decomposition by back 
substitution. 
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The authors are not familiar with the capabilities of the current generation of STAR 
machines or plans for the next generation. It may be that the increased capability of 
these machines could solve the problem as posed. Also, the capabilities of sparse matrix 
routines have not been thoroughly investigated. It would appear feasible, however, to 
actually obtain the inverse of [A(M,o>)] for two-dimensional problems of practical size 
but not for three-dimensional problems. 

The alternative would be to calculate the transonic air forces in terms of a limited set of 
reference modes. Then as the mass and/or stiffness distributions of the basic structure 
are changed, the new natural modes of the modified system would be found as a 
superposition of the reference modes. In the same fashion, the generalized air forces for 
the modified system would be obtained as a superposition of the air forces for the 
reference modes. These reference modes would usually be the natural modes of the basic 
configuration. 

The generalized forces for the reference modes (or any other set of modes), if the inverse 
of [A] is not available, would be calculated from a sequence starting with the solution 
for {<p\\ equation (21) using LU decomposition. It is noted that LU decomposition is 
done once and the results are stored so that obtaining! <£i}for new modes is relatively 
efficient. The column of generalized forces is then found from the matrix triple product 

rxl rxn nxn nx4n 4nxl 

jo s ) =m T ([T,] ITp]){ 

The{ is a subset of the full {<p \} matrix from the solution of equation (21). The 

matrix product enclosed in the parentheses may be calculated once and stored for future 
use. The integer r is equal to the number of reference modes used, which would be 
expected to be approximately 20 or 30. Generally, it appears that obtaining the <pi’s for 
additional reference modes will take less time than was required to do the LU 
decomposition. 
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11.0 CONCLUSIONS 


This report has further explored a particular finite difference formulation for analyzing 
unsteady transonic flow over harmonically oscillating wings. The preferred numerical 
solution process of line relaxation has proved unstable for ranges of Mach number (both 
subsonic and supersonic) and reduced frequency of direct interest in flutter analyses. 
Although no means were found for extending the range using relaxation procedures, a 
direct solution of the complete finite difference mesh was shown to produce solutions at 
subsonic Mach numbers outside the range of solution convergence for the relaxation 
process. It is surmised that the direct solution could also be applied to flows with a 
supersonic freestream Mach number as well. 

Because of limited computer capacity, the direct solutions obtained were for a coarse 
mesh, and accuracy was observed to decrease with frequency. The means for improving 
the accuracy of the direct solution are indicated by a study of a similar one-dimensional 
problem for which exact, analytic solutions were readily obtainable for comparison with 
the solutions from the finite difference analysis. The accuracy of the finite difference 
procedure was found to be proportional to n Xj so that the mesh spacing must be varied 
inversely to the 3/2 power of frequency if accuracy is to be retained. For the higher 
values of reduced frequency at values of Mach number close to 1, this will mean 
working with very large sets of finite difference points. How the use of nonuniform 
mesh spacing will affect this conclusion has not been examined. 

In addition to the general error level, large excursions in error are caused by the 
presence of real eigenvalues associated with the mesh region and the far-field boundary 
conditions. These excursions can be supressed in the one-dimensional examples by 
proper selection of boundary conditions that result in the replacement of the real 
eigenvalues by complex eigenvalues. Since these boundary conditions are in the nature 
of outgoing waves, it is assumed this can be done in the two- and three-dimensional 
analyses also. 

The three-dimensional program developed in reference 5 was extended to analyze wings 
with swept leading and trailing edges, and solutions for both a moderate-aspect-ratio 
rectangular wing and a low-aspect-ratio, clipped delta wing are presented. The pressure 
distributions appear reasonable although, as yet, no experimental results are available 
for correlation. 

A conservative shock point operator has been derived for use in two- and 
three-dimensional analyses. In the rectangular wing analyses, use of this operator 
significantly increases the effect of the shock on the unsteady pressure distributions and 
smooths the spanwise distribution of pressure. In relaxation calculations for a 
supersonic freestream, use of the shock point operator extends the range of convergence 
but does not remove the relaxation instabilities. 

Boeing Commercial Airplane Company 
P.O. Box 3707 

Seattle, Washington 98124 
November 1977 
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APPENDIX A 


ONE-DIMENSIONAL PROBLEM 

A.l PROBLEM STATEMENT 


To gain insight into the accuracy obtainable by using direct solution methods as well as 
the relaxation solution stability problem, a one-dimensional version of the flat-plate 
small perturbation equation was investigated. This equation, which is obtained simply 
by dropping the <plyy term from the two-dimensional equation, is 

v .. 2ico co 2 > 

— I/Jj +— y?i = 0, cj > 0 

i xx e e 1 


2 2 

or, dividing by K and using eK = (1 - M )/M 

<Pl - 2iAjMi£>i^ + Xj 2 (1 - M 2 ) = 0 (A-l) 

2 

where Xj = coM/(l - M ). The problem then was the solution of equation (A-l) for v?i(x) 
on an interval x = a to x = b with specified types of boundary conditions; i.e., a 
two-point boundary value problem. 

The problem was numerically solved by discretizing the derivatives with second-order 
approximations on a uniform mesh, as the two-and three-dimensional cases in reference 
5 were solved. The numerical problem thus becomes one of solving a linear system of 
the form = B, where A is a tridiagonal matrix and B depends on the boundary 

conditions. The solution was obtained with the tridiagonal solver used for each row of 
the two-dimensional row relaxation solution. 


A.2 ANALYTICAL RESULTS 

The advantage of experimentation on such a simplified problem is that analytic results 
are readily available for comparison with numerical calculations to obtain accuracy 
information. 

To begin, we observe that the general solution of equation (A-l) is given by 


^l(x) = 


c i 


iXi(l+M)x 

e 


+ C2 e 


-iXj(l-M)x 


(A-2) 


where C x and C 2 are independent of x and are to be determined by the boundary 
conditions. We note that in general <pi contains components of substantially different 
wavelengths, in fact in the ratio (1 + M)/(l - M) = 19/1 " for M = 0.9. Since 
approximation of the shorter wavelength component is less accurate for a given number 
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of mesh points than the longer wavelength component, we would expect the error to be 
larger when Cj is such that this component is significant, regardless of the value of X*. 
This was indeed found. (See fig. 4.) Since Cj and Cg are determined by the boundary 
conditions, we turn next to consideration of these. In their most general form, these are 

aj = 7j at x = a 

x (A-3) 

(*2 + 1 ® 2 ^ 1 = ^2 x = b 

where the P\, and yi are independent of x. Special cases of these conditions include 
Dirichlet (aj = a 2 = 0, Pi = p 2 =1), Neumann (a 1 = a 2 = 1, Pi =/3 2 = 0),and the third 
kind (ai = a 2 = 1, Pi ^ 0, p 2 ± 0). 

Since Cx and C 2 are determined from the boundary conditions, the choice of the 
constants in equation (A-3) will clearly affect the accuracy of the solution. There is a 
less obvious way in which the choice of these constants affects the accuracy: certain 
choices will lead to nonuniqueness of the solution, which is reflected as large increases 
in the error for certain reduced frequencies. To be more precise, for certain choices of 
a 1? Pi , « 2 , and p 2) there are values of Ax (eigenvalues) for which nonzero solutions 
(eigenfunctions) to equation (A-l) with yi = y 2 = 0 in equation (A-3) will exist. 
Information is easily obtained as to the values of the a and p constants in equation 
(A-3) for which real eigenvalues exist. Substitution of the general solution equation 
(A-2) into (A-3) with yi = y 2 = 0 yields a 2 by 2 linear system for Ci and C 2 . In order 
for this homogenous system to have a nontrivial solution (i.e., for eigenfunctions to 
exist), it is necessary that the coefficient matrix be singular (i.e., have a 0 determinant). 


When the substitution is made and the determinant set equal to 0, one obtains, after 
some simplification, the relation: 

[Aj (ajjf^-a^j)] cosXj (b - a) 

* | (1 - M~) + /3 1 ^2 1 + iAjM (aj^ + o^jSj) J s * n (b - a ) = 0 (A-4) 

Note that if either the coefficient of the cosine term or the sine term is 0, then real 
values of Xj exist for which equation (A-4) will be satisfied; for example, 

Aj ^ _ g > m ~ 1,2, (A-5) 


This will clearly be the case for Dirichlet {a\ = a 2 = 0, Pi — p 2 = 1) and Neumann (ax 
= ol 2 = 1, Pi — p 2 =0) boundary conditions, since the coefficient of the cosine term is 
then 0. The same is true for boundary conditions of the third kind (ax = a 2 =1) when 

Pi = P2-) 

There are other values of the a’s and p's which lead to real eigenvalues; for example, 
values such that the coefficients of both the sine and cosine terms are real or pure 
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imaginary. Of more interest, however, are values such that equation (A-4) cannot be 
satisfied by real Xj . Such values can easily be obtained using mixed boundary 
conditions. For example, specifying <p x at x = a and <pj X + ico(l - M)<p* at x = b implies 
that in equation (A-3) a 1 = 0, fi 1 = 1, a 2 = 1, and £2 = - M). Substituting these 

values into equation (A-4), we find that this equation becomes 

% iXi(b-a) 

- Xje 1 ^0 

which cannot be satisfied for real Xi ^ 0. Thus these boundary conditions do not permit 
real eigenvalues. Similarly specifying <pi x - iAjd + M)<pi at x = a and <pi at x = b 
implies <*i = 1, = -iXi(l + M), a 2 = 0, and /3 2 = 1, which when substituted into 

equation (A-4) yields the equation 

X,e iX ' <b - a) = 0 

which cannot be satisfied for real Xj ^0. Thus again the boundary conditions do not 
permit real eigenvalues. The difference in the behavior of the error for real and complex 
eigenvalues is illustrated in figure 5. 

Some idea of the error introduced by discretization of the analytical problem may, at 
least in the case of Dirichlet boundary conditions, be gained from the literature. In ref. 
18, Fischer and Usmani have shown that for the problem 

i//" + X ^ = o, xe [a, b] 

*Ka) - iKb )=^ 5 

when i//" is replaced by the usual second-order finite difference approximation on a 
uniform grid of spacing h, the maximum absolute error, E max , satisfied E max < Ex, 
where 

Ej = h 4 m 4 N/[ 12 sin 0 • | sin (N + 1) 0| ] 

where 


1U4 = 


max 
xe[a, b] 


(X) 


cos 6 = 1 h 2 \ j 2 


h = 


b-a 
N + 1 


When hXj is small, we have 0~hXi, and 



m 4 (b - a)/[ 12Xj 


%/l-^h 2 



IsinX, (b-a) |] 


Since tfj = D x sin Xi x+ D 2 cos Xix, we have nq = DXi 4 for some constant D, so that Ej 
behaves like Eh 2 Xj /| sin Xj(b - a)| , with E some constant. 
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o 

Now since equation (A-l) for <p i may be transformed into i/f" + Xj iff = 0 by the 
nonsingular transformation (p j it is not unreasonable to expect the error in 

the numerical solution of the discretized version of the (fi problem, and this in fact was 
found to be the case. (See fig. 6.) 

The results and conclusions from this analysis are presented in section 5.1. 
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APPENDIX B 


EVALUATION OF FAR-FIELD WAKE INTEGRAL 
FOR AN ARBITRARY WING 


We consider the evaluation of the field wake integral given in equation (B-l) of 
reference 5, namely, 


.+y t 


<pi (*i 


1 / ’ icoxt(yi') 

’ y i’ z i ) = ^T J e A ^i t (yi') d y T 


-y t 


/. 


Vvi') 

where the partial derivative is to be evaluated at 

-i\j [MCxj-xjO-R] 

\p = e 


-IWXi dxb . , , >u , 

e az^' (x i " x i 5 y i “ y i 5 z i “ z i )dx i (B-d 


and 

R=>/(xj -xj ') 2 + (yi -yi ') 2 + Z ] 2 

As before, the evaluation will be carried out for x 1 = 1.0. The trailing-edge function will 
be assumed to be a single-valued function of the form 

x te (yj') = 1 + f(y i ') 

where f(yi') >0 for -yt^yt'^Yt’ but * s otherwise arbitrary. For the rectangular wing, for 
example, f(yi') = 0, while for the swept wing with straight trailing edge, 
f(yj') = a - 1 y x ' | , where a is the tangent of the sweep angle. A more general example is 
given in Figure (B-l). 

Equation (B-l), after taking e 1Ct) ‘ 10 into the xi '-integral, becomes 


(i, y i 

1 w 


r +y t 

= z i )= 4 hJ 

-y t 


e 10jf(y| ) A^i te (yi')dyi' 


. f e -ico(x, '-1) 

X J '= l+f(y ! ') 


3z j 


(1 -xj', y] -yj',zi 


zj'ldxj' 


(B-2) 
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Let I w be the inner, x x '-integral 


A w 


/ oo 

e " iw(x i 

A,'=f(yi') 


_1) r z / (1 _x i'’ y i * y i'’ z i ” z i' )dXi ' 


:i'=f(yi') ^ 

setting p = Xj' - 1 and inserting the expression for rp yields 

-iXj (Rj+Mp) 


f °° * - iX l 1 

t _ / -iwp _<L .§ 

w 4,-> 3zi 1 


dp 


‘1 


with 


R 


=>/ p 2 + Rq 2 and R 0 2 = (yi -y 1') 2 + (zi -zi ') 2 


Taking the d/dzj 'outside and combining the exponentials, we obtain 


or since 






M(cj + X, M) 

f 00 e _] M" MR l + xT — P 

= JL_ / 5 1 dp 

3Z 1 - / f(y,') 


/ , coM 2 \ 

coM/(l - M 2 ) 


J CO + 

M(co + Xj M) _ 1 \ 1-M 


Next, let 


1 

W 3 Z J 


O 

/ 


-i^ ( M R i +p) 

e dp 


R 


fCy 1') 


U = 


MRj + p 

W 


When 


and as 


M^f 2 ( yi ') + R 0 2 + f(yr) 

P = f (Yl'X u= ^ = U 1 


p-*- °° u -*- 00 


(B-3) 


(B-4) 


(B-5) 


(B-6) 


(B-7) 


(B-8) 


(B-9) 
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Further 


/3 2 R 0 2 + M 2 (p 2 + R 0 2 ) + 2MRjp + p : 


n 9 

which becomes, using /3 = 1-M 


+ u 2 = ^ j^R 0 2 + P 2 + 2MRjp + M 2 p 2 J 


= ^[ R 1 2 + 2 MR 1 P + M2p 


2 1y 2 RjjHMp 


(B-ll) 


du _ 1 

dp |3R 0 


1 [ M p/R I + 1 ] =J-. 


1 R l +M P _ Vl+U 2 
Rj /3R 0 Rj 


(B-12) 


Xj0 R o 


,3 £_ — M_ 


w 0 Z j 


uj V 1 + u" 


Now, let u = y + uj so that 


Xj^Rq 

a _i ~ M~ (v+u 1 } H 

T = J)_ f e 1V1 Av 

w " 3z ll ./l+fp + u,'/ 


1 •'o v 1 + + 11 1) 

The singularities of the integrand are where 

1+ 0 + Ui) 2 = o 


(B-13) 


(B-14) 


(B-15) 


v = -u i ± i 


Now from equation (B-9) and the assumption that f(yi') > 0, it follows that uj>0, so the 
singularities are in the left half of the complex plane. Thus, applying Cauchy’s theorem 
to the contour integral 


J e M 1 dv 

yfi + O' + uj) 2 


(B-16) 
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where the contour is shown in the following sketch 



- Re (*0 


Then, c = cj + C2 + C3 and 


V c ! W C +I W C +I W C 0 
c Cj C2 C3 


and 


T = 

w r -+oo^- ' 


d T _ lim T 9_ T d . ~| 

9 z l' C 1 r^.ool_ & z i' l °2 3 Zl ' ic 3J 

As r-^oo, then integral I WC2 ~>0. On the contour C 2 , v = re 

-iX^R^uj/M) p/2 e -iXig(R 0 /M) r expC-jg^ g 

J 0 yr+ 


and 


! w =' ire 
c 2 


(r exp(-ifl) + u j)^ 

For r sufficiently large, the denominator of the integrand is >r/2. 


dd 


i w 1 <2 r /2 e -M« R o/M) r si„ 6de 

c 2 •'O 


/ 


<2 1 e 

,/ 0 


n/2 -(2/7r)X 1 (R o /M)-r0 


d0 


7rM 


Xl/3Ro 


-[,. e -MW-/M)]^ 0asr . 


(B-17) 


(B-18) 
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Thus, we have 


T =- ^ T 

r->oo3 Zl f W G 3 


A w 


(B-19) 


On Cz,v = irj and dv = id r\ so that 


or letting 17 = -77 


V =i 


Jf 


X]/3R 0 

0 _i Cir?+u 1 ) 


dr? 


w, 


c 3 J r \/l + (ii? + u j Y 


r r -i— M~ (u l' ir?) , 

= j j £ d?7 

•'a 


c 3 


L) >/ 1 + ( u i - ir ?) 2 


(B-20) 


(B-21) 


The next step is to find the partial derivative in equation (B-19). We note from equation 
(B-9) that u x is a function of Ro and therefore of zj ' and thus we must also differentiate 
Uj. It appears simplest to move the u j dependence to the limits before differentiating. 
Letting /x = Uj - irj, then d/u. = idr?, or d/x = -idrj, and when tj = 0 , /x = u^; when 17 = r, 


/x = Uj -lr, so, 


w 


c 3 


-/ 


Ui-ir -iXjj3Rg(p/M) 


x/l + 


dp 


dzj' w, 


_( /9u 1 ^f e _iX l^ u r ir) / M ^^RqUj/M- 

' S 1 L - 1 ^ 2 ' - 


/ 


u r ir 3 -iXj/3RQp/M 


9zi ' 


1 J\ +/x : 


dMi 


(B-22) 


Performing the differentiation, the integral becomes 

^ui-ir n ( i\K\ -iXjlSRoM/M 
f 1 -iAj/3 (/i/M) e 1 u 3 Rq 

J ^=5 


u 


\/l +M 2 

or, transforming back to tj 

r-iX,4 (u 1 -W/M]e- iX l SR 0 <U l- 1 ’ ,)/M 


dzj' 


- 1 




\/l + (u j - ir?) 2 


dr? 


3Rq 
3zi ' 
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or 


ft -iAi0 R O Ui/ M f T u l~ iT ? -XtfRtfi/M 3Rp 


M 


r u i 

i yr + 


(uj -ir)Y 


9zj ' 


(B-23) 


Substituting equation (B-23) into equation (B-22) and taking the limit as r we have 

d T _ -iXi^R nU ,/M[~ i /duA 

dz >' Wc 3 

Xj/3 f°° uj-iT? -X,)3R 0 r?/M dR 0 
+ m J , - e dr? • , 

0 v 1 + (u j - ir)Y 


(B-24) 


A 

Let I w be the 


then 


integral in equation (B-24). Integrating by parts with 

-X^R n r?/M . uj-ir? 

u = e 1 u dv = — 1 ■ ■ 


y/i + (uj - i t\Y 


du = - Xj|3(R 0 /M) e X 1^ R 0^/ M dr? 


and 


from which 


' = + (uj -i??) 2 


^^e-X^Ro./My^— -5 ^ + iMRo f* 

t?=o * / n 


e -X^R 0 T?/M ^ 


= -i>/l + u i 2 + i ~ J \A + ( u 1 - iT ?) 2 e 

0 

Making the transformation r = Xi/SRqtj/M, we have 

I w = -iv/l+u, 2 + i \X + ^U| = -r 


dr 


I =— 

A w p 


-0\A +uj 2 + y 

0 1 u 


e~ T dr 


(B-25) 


65 



The final calculation required is that of dd 1 ldz 1 ' . Since u 2 depends on z x ', only through 
Ro, duj/dzx' = (dux/dRo) ( dR 0 /dzx ')• 


From equation (B-9), 


from which 


U 1 = 


M s /\ " 2 (y j ') + R 0 2 + f(y i ') 

0 R O 



(B-26) 


Finally, substituting from equations (B-25) and (B-26) into equation (B-24), and using 
equation (B-19) and the fact that dRo/dzj' = (zi - zjVRo, we have for zj' = 0, 


I 


w 


■iX |/3 Rqu j/M 


z j e 

1 

( M 

R o 


^v^Cyp + Ro 2 


R o 


l h 

M 


-$J\ + u j 2 + f v4 2 + ^ uj - ^r^ 2 e r dr 


(B-27) 


where the bars signify evaluation at zx' = 0. That is, 


R, 


o=\/& 7 - 


yf ) 2 + z r 


and 


_ _Myf 2 (y,') + R 0 2 + f(y 1 ') 
= (*0 


The integral in equation (B-27) differs only slightly from that given in reference 5 and 
may be evaluated by the same method, Gauss-Laguerre integration. 


Thus evaluation of the far-field wake integral is reduced to using equation (B-27) in 
equation (B-2) with the appropriate expression for the trailing edge. 


Specialization of equation (B-27) for the straight swept trailing edge is immediate, 
consisting only of taking f(yj') = a|yi'| where a is the tangent of the sweep angle 
measured back from the y 1 / axis. 

A check is available by comparison with the rectangular wing case. Taking f(yi') = 0, 
we have 


U| - M/P and 



UP 
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so that equation (B-27) becomes 


T = iL,- iX i R ol^i 

R 0 e ( M 

as was previously obtained. 



Acknowledgement. Calculations performed by R.W. Call toward extending the 
treatment in reference 5 to the case of the sweptwing with straight trailing edge have 
provided useful insights for the general treatment given here. 
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APPENDIX C 


DERIVATION OF SHOCK POINT OPERATOR 


When a rapidly decelerating flow is supersonic upstream of the point i, j (ui.i£,j<0) and 
the flow becomes subsonic downstream of the point Ui +1 /2,j > 0)> a shock wave then lies 
close to the point i, j. To satisfy the appropriate jump conditions across the shock, the 
difference operator for the point i, j must be conservative. To obtain such an operator, 
we apply the divergence theorem 

/ V-Fdv= f F.f?d s (C-l) 

*v J S 

to the differential operator expressed in conservation form for the control volume 
consisting of lines drawn midway between consecutive columns and rows of mesh points 
as shown in figure C-l. Here n is the outward normal to the closed surface. We shall 
consider only two-dimensional flow, but the generalization to three-dimensional flow 
requires only the addition of the central difference operator for <pi zz at the point i, j, k. 
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Figure C-1.— Control Volume for Shock Point Operator 
The basic partial differential equation, equation (17) of Ehlers (ref. 1), has the form 

— 2icos^j /e^ = 0 

Hence, the vector F is 

F = u^j /e, 
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and the approximate evaluation of the surface integral in equation (1) yields 




“l.*x , ~ u . 3^x o ~ 2kj /V I."*’. 3.V 6 

1+ 2 J i+J 1 I i~|j V 1+ 2 J 1 2y J 

+ j - ^x^y j + ^ij ^x^y^ij “ ^ 


where i + M: denotes the value of the quantity at the point midway between Xj and xj +1 ; 
and similarly, for the other half integer subscripts. Dividing the equation by 9.Jt y puts 
the equation in difference form: 


1 " U 3.^X n 

1+ 2 J i+^j '~2 j i-|j_ 


/2 X - 2(iw/e) 


^.,1.-^. 3. 
1+ 2 J 1_ 2 J 


/«> 


y... i i 
U+2 


/Cy + Qij^ij - 0 


(C-2) 


Now 


V. 3. V C x ~ \ X.-*. 1 . + 1 . " 3 . ) / C 




1+3J i-fJ 1_ 2J 1_ 2 j / 


( x i+l - x i-l ^Xy + ( x i - x i-2^ 


/ x l- x 3 
\ 1 2 1 2 / 


x i+l + x i - x i-l “ x i-2 


We substitute the central difference operator for <p Xm . and upwind difference operator for 

<£*.. and obtain IJ 

A ij 

= E)Xj jc 2 j(v ? ij - ^i-lj) " d ii-1 (v?i_ij -^i_2j|j/( DX l + DX 2 ) 

+ DX 2 [c h (y> i+lj - - y?ij) + d u Opy - /(DX t + DX 2 ) (C-3) 

where DXj = xj - X[_ 2 , DX 2 = Xj +1 - x^-*, and c n, d^, and c 2 j are given in equations 
(19), (20), and (26) of reference 1. 

( u .A.*x , - u . 3.^x o V*x 
\ 1+ 2 J i+ 2 J 1 2 J i~|j / 

= D 1 f* 2c M u j . (<Py - - 2d M u 3 - ^_ 2j )l 

L i 2 J 1 2 J 

+ D 2 r 2Ci u (*i +lj - <Py) - 2dj U _J (spy - (C-4) 

L i ' 1 oJ J 
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with D x = DXi/CDXi + DX 2 ) and D 2 x DX 2 /(DXx + DX 2 ). The y derivative becomes 

*Vyy = (*V 1 " *V„ i V £ y = 2 a j^-i - 2 < a j + b j> m + 2b j ^j+i (C . 5) 

\ y+2 y 2 ' 

where aj and bj are given in equation (23) of reference 1. 

Finally the difference equation at the shock point ij is obtained by substituting 
equations (C-3), and (C-4), and (C-5) into equation (C-2). After some simplification, we 
obtain 


ajY>jj_i - (aj + bj + D 2 Ej + D 2 E 2 - Dj E 3 - qy/2) spy + bj ^ ij+] - 

(DiE 3 + DjE 4 - D 2 E 2 ) - DjE 4 D 2 E 1 y’i+lj (C-6) 

where 

E) = c i u j - icoc yje 

i+ 2>. 

E 2 = dj u j +icodjj/e 

i- 2 j 

E3 = Cj_2 11 1 - icoc2j/e 

i- 2 j 

E 4 = d i-1 u . 3.- icod li-l l £ 

X -2 J 

When Dj = 1 and D 2 = 0, the difference equation (C-6) reduces to the hyperbolic 
upwind difference equation; and when D 2 = 1 and Di = 0, it reduces to the conventional 
elliptic central difference equation. For equally spaced points, Di = D 2 = %. With this 
value for T>i and D 2 , equation (C-6) becomes the shock point operator used by Traci, 
Albano, and Farr (ref. 9). When Dj = D 2 = 1, we obtain the shock point operator as 
used by Murman (ref. 23). The present formulation as well as that of reference 9 yields 
a form of the differential equation consistent with that at adjacent points, whereas 
Murman’s does not. 
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APPENDIX D 


RESIDUAL ANALYSES 


This appendix derives an interpretation of RESIDUAL as a mass flux and a 
mathematical explanation of the reason why the RESIDUAL and ERROR curves as 
plotted in a convergence history are essentially parallel. It also includes numerical 
examples illustrating these concepts. The differential equation that we solve for the 
unsteady flow potential is basically the equation for the conservation of mass at a point 
in the flow. This is easily seen since it results from separating the continuity equation 

P t /P + V m Vyp\ +V<Pi *Vp/p = 0 

into steady flow and unsteady flow equations after Vp/p is eliminated by using 
Bernoulli’s equation (see page 35 of Ehlers, ref. 1). The differential equation for the 
unsteady complex potential can be expressed in the form (eq. (75) of ref. 1): 

V* F + = 0 

where Q does not depend upon the unsteady potential <pi, and the factor F= (F 1 ,F 2 ,F 3 >. 

Consider now a rectangular parallelepiped of sides Ax, Ay, Az centered about the point 
x,y,z. The mass produced inside the region is found by integrating the differential 
equation over the volume. Applying the divergence theorem to the first term of equation 
(19) gives us 

AyAz [Fj(x + Ax/2, y, z) - F ] (x - Ax/2, y, z)] 

+ AxAz [F 2 (x, y + Ay/2, z) - F 2 (x, y - Ay/2, z)] 

+ Ay Ax [F 3 (x, y, z + Az/2) - F 3 (x, y, z - Az/2)] 

+ AxAyAz Q(x, y, z) i/?j =0 
Factoring out AxAyAz yields 

AxAyAz | [ F j (x + Ax/2, y, z) - F j (x - Ax/2, y, z)] /Ax 

+ [F 2 (x, y + Ay/2, z) - F 2 (x, y - Ay/2, z)] /Ay 
+ [F 3 (x, y, z + Az/2) - F 3 (x, y, z - Az/2] /Az 
+ Qsf 5 1 1 = 0 

The quantity in brackets can be recognized as the difference equation which we are 
solving for the unsteady potential. The complete left-hand side of the previous equation, 
which may not be necessarily zero, is the residual as used in this report and is a 
measure of the excess (or deficit) of the mass flux within each mesh in the flow field. 
Since for an exact solution of the difference equation this residual should be zero, it is a 
measure of how close the relaxation solution is to being converged. 
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A mathematical explanation of the reason why ERROR and RESIDUAL curves should 
be essentially parallel is also readily available. Block (i.e., row or column) successive 
overrelaxation for the solution of A <p{ = R may be written in matrix notation as 

^ 1 ( m+ D= (D + rL) -1 (-rU + ( 1 - r)D)i^ m ) + r ( D + rL ) _1 R 

where D consists of the blocks on the diagonal of A; L and U consist of the blocks below 
and above the diagonal, respectively; and r is the relaxation factor. (The special case 
where A is Hermitian was discussed in ref. 5.) Replacement of U by A-D-L followed by 
some matrix algebra yields the equivalent form 

^ ] (m+l) = ^ i (m) + r(] + rD -l L) _1 D -1 (R_ 

Here the last factor on the right is immediately recognizable as the vector of (unsealed) 
residuals. Multiplication of the residual vector by D" 1 effectively scales the residuals by 
an area, similar to the scaling used in the program, since the primary components of 
the elements of D have the dimensions of (area)” 1 . Moving the first term on the 
right-hand side to the left, dividing by r, and taking the maximum norm of each side, 
we have 

~ll ^( m+ l ) - v?j(m) || <|| (H-rD^L)" 1 || • || D -1 (R-A^/ m >)|| 

Now, the left-hand side is just ERROR, as previously defined, while the right-hand side 
is just a constant times the maximum scaled residual. When this inequality is viewed as 
providing an estimate of ERROR in terms of the scaled residual, it is clear that these 
two quantities are roughly proportional. The essentially parallel nature of the ERROR 
and RESIDUAL error curves on a semi-log plot in both convergent and divergent cases 
is thus natural and expected. 

Examples of the behavior of the maximum RESIDUAL with iteration for a 
two-dimensional configuration are given in figures D-1 through D-7. Two different 
maximum RESIDUALS are shown. The first, called an "INTERMEDIATE RESIDUAL” 
is calculated after each line relaxation is completed. The INTERMEDIATE RESIDUAL 
is computed using both old and new values of the velocity potential and, if the 
relaxation factor is set to 1.0, it should be zero. In our calculations it is of 
approximately 10” 15 for an ORF of 1. It is of interest mainly because it can be 
calculated with the coefficients and velocity potentials available at the end of each line 
relaxation calculation and thus may be obtained very efficiently. The second line shows 
the true maximum RESIDUAL, which is calculated in a separate pass after all velocity 
potentials have been updated. It includes the relaxation factor and thus reflects how 
well the calculated velocity potential as scaled with the relaxation factor satisfies the 
finite difference equation. 

In figure D-1, results are presented for a flat plate example (no mixed flow) solved with 
row relaxation using an ORF of 1.85. As would be expected from the preceding 
discussion, the true RESIDUAL curve runs parallel to the ERROR curve, and in this 
particular case the true RESIDUAL curve lies on top of the ERROR curve. The 
INTERMEDIATE RESIDUAL curve lies above the other two. Figure D-2 presents the 
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Figure D-3.— Sample ERROR and RESIDUAL Curves Versus Iteration for Row Relaxation 
ORF * 1.85 , M = 0.9, c o = 0. 12, Flat Plate 


i 












3.3E-6< ERROR < 6.6E-6 
1.0E-6 < ER ROR < 3.3 E-6 
6.6E-7 < ERROR< 1.0E-6 
3.3E-7 <ERROR < 6.6E-7 
ERROR <3.3E-7 


Figure D-5.— Areas of Matching ERROR Values for a Converged Solution, Row Relaxation, 
ORF = 1.85, M = 0.9, co = 0.06 
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3.3 E-6 < RESIDUAL< 6.6E-6 
1.0E-6 < RESIDUAL <3.3E-6 
6.6E-7 < RESIDUAL < 1.0E-6 
3.3E-7< RESIDUAL < 6.6E-7 
RESIDUAL <3.3E-7 


Figure D-6.— Areas of Matching RESIDUAL Values for a Converged Solution, Row Relaxation, 
ORF = 1.85, M = 0.9, co = 0.06 
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Downstream boundary 


Upstream boundary 



Figure D-7.— Areas of Matching ERROR Values fora Converged Solution Row Relaxation, 
ORF = 1. 70, M = 0.9, w = 0.06 
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Downstream boundary 


results for the same solution but with ORF = 1.70. This time there is a distinct knee in 
the curves at about 50 iterations. After the knee in the curve, the true RESIDUAL 

lies between the INTERMEDIATE RESIDUAL and ERROR. For this particular case 
and with the level of convergence set at 10‘ 5 , an ORF of 1.7 would be better than an 
ORF of 1.85. However, if the convergence criterion were set lower, then the ORF of 1.85 
would be preferred. 

For figure D-3 the problem has been changed by raising the frequency so that the 
solution diverges. Again all three curves move together but in a far from smooth 
manner. Both the ERROR and the INTERMEDIATE RESIDUAL provide a good 
indication of what is happening to the solution. Figure D-4 presents ERROR and the 
true RESIDUAL for a case with thickness, and again the two curves are essentially 
parallel to each other. 

Figure D-5 shows the ERROR distribution for a converged (maximum error less than 
Iff 5 ) solution. The finite difference area below the wing is divided into levels of 
ERROR. For example, the most converged points in the field according to the ERROR 
criterion are in the uniformly shaded areas adjacent to the wing leading edge and near 
the lower boundary beneath the wing. Figure D-6 shows the same type of plot for the 
RESIDUAL. Here, the most converged portion of solution according to the RESIDUAL 
criterion lies nearly uniformly under the wing extending from wing surface to lower 
boundary. Figure D-7 shows the ERROR distribution for a different overrelaxation 
factor (1.7 rather than 1.85) with the distribution quite different from that of figure 
D-5. 
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